arXiv:1509.05232vl [cond-mat.stat-mech] 17 Sep 2015 


The critical catastrophe revisited 


Clelia de Mulatier 

CNRS - Universite Paris-Sud, LPTMS, UMR8626, 91405 Orsay Cedex, France 
CEA/Saclay, DEN/DANS/DM2S/SERMA/LTSD, 91191 Gif-sur-Yvette, France 

Eric Dumonteil 

CEA/Saclay, DEN/DANS/DM2S/SERMA/LTSD, 91191 Gif-sur-Yvette, France 

Alberto Rosso 

CNRS - Universite Paris-Sud, LPTMS, UMR8626, 91405 Orsay Cedex, France 

Andrea Zoia 

CEA/Saclay, DEN/DANS/DM2S/SERMA/LTSD, 91191 Gif-sur-Yvette, France 
E-mail: andrea.zoiaOcea.fr 

June 2015 


Abstract. The neutron population in a prototype model of nuclear reactor can 
be described in terms of a collection of particles confined in a box and undergoing 
three key random mechanisms: diffusion, reproduction due to fissions, and death 
due to absorption events. When the reactor is operated at the critical point, and 
fissions are exactly compensated by absorptions, the whole neutron population 
might in principle go to extinction because of the wild fluctuations induced by 
births and deaths. This phenomenon, which has been named critical catastrophe, 
is nonetheless never observed in practice: feedback mechanisms acting on the 
total population, such as human intervention, have a stabilizing effect. In this 
work, we revisit the critical catastrophe by investigating the spatial behaviour of 
the fluctuations in a confined geometry. When the system is free to evolve, the 
neutrons may display a wild patchiness (clustering). On the contrary, imposing a 
population control on the total population acts also against the local fluctuations, 
and may thus inhibit the spatial clustering. The effectiveness of population control 
in quenching spatial fluctuations will be shown to depend on the competition 
between the mixing time of the neutrons (i.e., the average time taken for a particle 
to explore the finite viable space) and the extinction time. 
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Many physical and biological systems can be represented in terms of a collection of 
individuals governed by the competition of the two basic random mechanisms of birth 
and death. Examples are widespread and encompass neutron multiplication [Dials], 
nuclear collision cascades ISIEIEI, epidemics and ecology la H n 0 HU], bacterial 
growth HD Ha, and genetics Ha HI HSj- Neglecting particle-particle correlations 
and non-linear effects, the evolution of such systems can be effectively explained by 
the Galton-Watson model [3]. When the death rate is larger than the birth rate, 
the system is said to be sub-critical: the population size decreases on average, and 
the ultimate fate is extinction. This occurs for instance for nuclear collision cascades, 
where charged particles are progressively scattered and absorbed by the medium [3 [5]. 
When on the contrary the birth rate is larger than the death rate, such as for bacteria 
reproducing on a Petri dish [12], the system is said to be super-critical. In this case, 
the population size grows on average. However, because of fluctuations on the number 
of individuals in the population, a non-trivial finite extinction probability exists for the 
whole system |3] . A super-critical regime is typically found also during the early stages 
of an epidemic (the so-called ‘outbreak’ phase), where a fast growth of the infected 
population is observed, until non-linear effects due to the depletion of susceptible 
individuals ultimately slow down the epidemic [T6] . 

In the intermediate regime, the population stays constant on average, and the 
system is said to be exactly critical. A prominent example of a system operating 
at the critical point is provided by the self-sustaining chains of neutrons in nuclear 
reactors HI 0 HD- Avalanches and self-organized criticality are other examples of 
system operating at or close to the critical point. At the critical or nearly critical 
regime, fluctuations due to birth and death may become particularly strong [3]- 
Assuming without loss of generality that the birth and death rates are Poissonian, 
from the Galton-Watson theory it is known that at criticality the total number N(t) 
of particles in the system stays constant on average, i.e., {N) = Nq, whereas the 
variance grows in time, i.e., Var[N] = {N‘^) — {N)‘^ oc XNot, where A is the birth/death 
rate |3|- This immediately implies that the typical fluctuations of the population 
size, say ^ Var[N ], will become comparable to the average population size Nq over 
a time ^ Nq/X. Hence, a critical system will have a characteristic extinction time 
of the order of te — Nq/X This quite extreme behaviour has been observed for 
instance in avalanche dynamics. In the context of neutron multiplication, the shut¬ 
down of a reactor operated in the critical regime due to the extinction of the fission 
chains has been theoretically investigated [DEHaHlI and goes under the name of 
critical catastrophe [2]. However, such observation is in open contradiction with the 
experience: the behavior of nuclear reactors at the critical point is actually stable. 
This apparent paradox has been explained by pointing out that including feedback 
mechanisms (representing for instance human intervention) in the Galton-Watson 
model induces a stabilizing effect acting against the total population fluctuations [2|- 

In most of the models cited above, individuals also interact with the surrounding 
environment and are typically subject to random displacements millinillQ]. The 
interplay between the fluctuations stemming from birth-death events and those 
stemming from diffusion will thus subtly affect the spatial distribution of the particles 
in such systems |20l Ell EH ESI El] • In particular, it has been shown that at and close 
to the critical point a collection of such individuals, although spatially uniform at 
the initial time, may eventually display a wild patchiness (see Fig. [^, with particles 
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Figure 1. Monte Carlo simulation of the evolution of a collection of particles 
in a two-dimensional box. The particles are prepared at t = 0 on a uniform 
spatial distribution. In case a), particles follow regular Brownian motions: as time 
increases, positions are shuffled by diffusion, but the spatial distribution of the 
particles stays uniform. In case 6), particles follow a branching Brownian motion 
with equal birth and death rates: as time increases, the population undergoes 
large fluctuations, and the particle density displays a wild patchiness. Eventually, 
the entire population goes to extinction. 

closely packed together and empty spaces nearby 13 ini Eg E]. Spatial clustering 
phenomena have been first identified in connection with mathematical models of 
ecological communities m HZ], and since then have been thoroughly investigated 
for both infinite and finite collections of individuals [9l [THl |25l |28l |29l |30j. Non- 
uniform neutron densities in the reactor might lead to local peaks in the deposited 
energy (hot spots) and represent thus a most unwanted event with respect to the safe 
operation of nuclear power plants [ST]. In view of the findings concerning the impact 
of a feedback on the global neutron population, we might wonder whether imposing a 
global feedback on the total neutron population affects also the local spatial behaviour 
of the particles. 

In this work, we will revisit the critical catastrophe of neutron chains in a 
prototype model of a nuclear reactor, with special emphasis on the spatial distribution 
of neutrons in confined geometries. We will show that actually a global feedback has 
a stabilizing effect also on the local particle density, and may thus inhibit clustering. 
In particular, the effectiveness of population control in quenching spatial fluctuations 
will be shown to depend on the competition between the mixing time of the neutrons 
(i.e., the average time taken for a particle to explore the finite viable space) and 
the extinction time. This paper is organized as follows: in Sec. we will sketch the 
statistical model of the neutron population in a nuclear reactor. Then, in Sec. 
we will illustrate the main findings concerning the pair correlation functions for a 
collection of particles in a confined box with and without population control. The 
detailed derivation of the pair correlation functions will be discussed in Sec. and 
conclusions will be finally drawn in Sec. Calculations for one-dimensional systems 
will be worked out in the Appendix. 
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Figure 2. Simplified scheme of neutron propagation within a nuclear reactor. A 
fission chain begins with a source neutron (marked with S in the figure) born from 
a fission event in the fuel. The neutron diffuses in the water and may eventually 
come back to a fuel element. Then, it can either be absorbed (these events are 
marked by magenta circles), in which case the trajectory terminates; or it can give 
rise to a new fission, upon which additional neutrons are set free (these events 
are marked in green), and the fission chain is kept alive. The system is operated 
at the critical point when the average net number of neutrons produced at fission 
is exactly compensated by the losses by absorptions. A control rod can absorb 
the excess neutrons so as to adjust the total population and enforce the critical 
regime. 


2. A prototype model of nuclear reactor 

A nuclear reactor is a device conceived to extract energy from the fission chains induced 
by neutrons [32]. To fix the ideas, here we will focus on the widely used light-water 
reactors. The nuclear fuel is made of uranium, arranged in a regular lattice and 
plunged in light water. A fission chain begins with a neutron emitted at high energy 
from a fission event (see Fig.[^. The neutron enters the surrounding water, slows down 
towards thermal equilibrium, and then starts diffusing. If the neutron eventually re¬ 
enters the fuel, it may i) be absorbed on the isotope of uranium, in which case 
the chain is terminated; or ii) give rise to a new fission event by colliding with the 
fissile isotope, whereupon a random number of high-energy neutrons are emitted. 
The water surrounding the fuel lattice acts as a reflector and prevents the neutrons 
from escaping from the core. A number of control rods are inserted into the core, 
with the aim of absorbing the excess neutrons and keep the population constant (this 
ensures a constant power output). When the neutron population grows, the control 
rods are inserted more deeply into the core, slowing down the chain reaction. On the 
contrary, when the population decreases, the control rods are raised, accelerating the 
chain reaction. 

Nuclear reactors are complex devices: their energy- and spatial-dependent 
behaviour can be fully assessed only by resorting to large-scale numerical simulations 
including a realistic description of the heterogeneous geometry m- However, for the 
purposes of this work we will introduce a simplified prototype model of a nuclear 
reactor that yet retains all the key ingredients of a real system. 
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We will assume that the reactor can be represented as a collection of N 
neutrons undergoing diffusion, reproduction and absorption within an homogeneous 
d-dimensional box of finite volume V, with reflecting (mass-preserving) boundaries. 
It is reasonable to require that the initial neutron population has a uniform spatial 
distribution. The stochastic paths of neutrons are known to follow position- and 
velocity-dependent exponential flights missiEii- For our model, we approximate these 
paths by regular d-dimensional Brownian diffusion with a constant diffusion coefficient 
D. The diffusing walker undergoes a birth-death event at rate A: the neutron 
disappears and is replaced by a random number m of descendants (m = 0 accounting 
for absorption), distributed according to a law Pm with average In 

order for the reactor to be exactly critical, we must have ui = 1. 

Clustering phenomena have been mostly analysed either in the thermodynamic 
limit {V ^ oo and N oc, with finite N/V [26j[29]) or in unbounded domains with 
finite N laiTo]. A realistic description of actual physical systems demands however 
that the effects due to the finiteness of the viable volume V be explicitly taken into 
account. In a previous work, we have shown that a neutron population within a finite- 
size reactor at the critical regime will ultimately undergo wild spatial fluctuations over 
the entire volume [30]. In this paper, we include the effects of population control by 
imposing that the total number N of neutrons in V is preserved, and investigate the 
consequences of such constraint on spatial fluctuations. The simplest way to enforce a 
constant N is to correlate reproduction and absorption events [3 HU]: at each fission, 
a neutron disappears and is replaced by a random number m > 1 of descendants, 
and m — 1 other neutrons are simultaneously removed from the collection in order to 
ensure the conservation of total population (see Fig. [^. This mechanism has been 
first introduced in the theoretical ecology (with binary branching = Sm ,2 OllOl), 
where similar large-scale constraints such as limited food resources have been shown 
to quench the wild fluctuations in the number of individuals that are expected for an 
unconstrained community. Similar effects have been also considered in the context of 
cellular growth under the effects of chemotaxis [35] . 

3. Physical observables and main findings 

Let us denote by n(xi,t) the instantaneous density of neutrons located at at time 
t. For an exactly critical system (z^i = 1), the average neutron density at a point x^ 
reads 

= Np{xi,t), (1) 

where we have set 

/9(Xi,t) = y dxoQ(xo)0(xi,xo,t). (2) 

Here Q is the spatial probability distribution function of the neutrons at time t = 0, 
and the propagator ^(x, xq, t) times dx gives the probability to find a regular Brownian 
particle inside the volume element (x; x + dx) at time t, knowing that the particle was 
at Xq at time t = 0. The propagator satisfies the diffusion equation 

^ 0 ,t)= (x, Xo, t), (3) 

with the appropriate boundary conditions. Note that at criticality the average neutron 
density becomes indistinguishable from that of N regular Brownian particles. This 
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Figure 3. The evolution of a collection of branching random walks with binary 
fission and = 1. At time t = 0, A = 3 particles are present, and the system is 
observed at successive times ti < t 2 < ts < t 4 . Top. When population control is 
not enforced, the number of particles present at the observation times fluctuates 
because births and deaths occur at random instants. Bottom. When population 
control is enforced, at each fission event a neutron is randomly chosen and 
removed, which exactly preserves N: at observation times, the total population 
is constant. 

result holds true independently of whether population control is applied. For an initial 
uniform source of particles Q = l/V^ it immediately follows that for a collection of 
N critical branching Brownian motions we simply have a uniform average density 
(n(xi,t)) = N/V^ at any time. 

In order to probe the spatial inhomogeneities of the neutron population due to 
clustering, we must therefore go beyond the average behaviour. A fundamental tool is 
provided by the two-point (or pair) correlation function h(x^,xj,t) between positions 
x^ and Xj, namely, the average density of pairs with the former particle in x^ and the 
latter in Xj. This quantity is proportional to the joint probability density for x^ and 
Xi m- For N independent random walkers (in absence of branching and death) we 
simply have 

= N{N - l)p(xi,t)p{xj,t). (4) 

In particular, if the particles are uniformly distributed, = N{N — l)jV‘^. 

More generally, the spatial shape of h((Ki^^j^t) conveys information on the correlation 
range, whereas its amplitude is proportional to the correlation strength. A flat shape 
implies that the correlations have the same intensity everywhere; on the contrary, the 
presence of a peak at x^ Xj reflects the increased probability of finding particles 
lying at short distances, which is the signature of spatial clustering igilolllSlEgllll]. 
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A closely related quantity is the average square distance between particles, i.e., 

2w^ ^ / dxi f dxj I Xj - Xjf Xj , t) 

J d-Kidxih{-Xi,:x.j,t) 

which is to be compared to the ideal average square distance of an uncorrelated 
population uniformly distributed in the available volume, namely, 

('^^hd j j 

where d denotes the spatial dimension. Deviations of {r‘^){t) from the reference value 
allow detecting spatial effects due to clustering j9l[TO]. 

Analysis of the model detailed above shows that the population dynamics is 
governed by two distinct time scales: a mixing time td oc jD and an extinction 
time te oc NjX. The quantity te physically represents the time over which a particle 
has explored the finite viable volume V by diffusion. Observe that the emergence of 
the time scale te is a distinct feature of confined geometries having a finite spatial size: 
for unbounded domains, te ^ oo. The quantity te has a different meaning according 
to whether population control is imposed lann]. For a free system, te represents the 
time over which the fluctuations due to births and deaths lead to the extinction of the 
whole population. For a constrained system with constant N, te represents the time 
over which the system has undergone a population renewal, and all the individuals 
descend from a single common ancestor. When the concentration N/V of individuals 
in the population is large (and the system is spatially bounded), it is reasonable to 
assume that te > te- Intuitively, the precise shape of the pair correlation function 
must then depend on the subtle interplay of te and te, where the former conveys 
information on the space exploration and the latter on the reproduction mechanism. 
Moreover, the pair correlation function will depend on whether population control is 
applied or not. In the following, we will denote by h/(x^,Xj,t) the pair correlation 
function for the case without population control, and hc(x^,Xj,t) for the case with 
population control. 

For the simple reactor model detailed above, the pair correlation function 
h/(x^,Xj,t) for the free critical system has been derived in a previous work [30] and 
is briefly recalled in Eq. in Section In this work we have explicitly derived the 
expression of hc(x^, Xj, t) for the case of a constant population in a confined geometry: 
the resulting formula is given in Eq. and the derivation is provided in Section 
Once Xj, t) has been determined, it is customary to introduce the normalized 

and centered pair correlation function, in the form 






(7) 


{n(xi,t)){n(xj,t)) 

which allows more easily comparing the amplitude of the typical spatial fluctuations 
to the average particle density. 

Eor the purpose of physical analysis and illustration, let us consider an initial 
collection of A" = 10^ neutrons subject to diffusion, reproduction and death in a one¬ 
dimensional bounded box [—L, L], with L = 1 (thus, V = 2). To fix the ideas, we will 
set a diffusion coefficient D = 10“^ and a birth-death rate A = 1. Eor this system, 
the mixing time reads te — 40.5 and the extinction time reads te — 10^ (hence 
see Appendix Appendix A). The initial condition for the neutron population 
is a uniform spatial distribution on [—L, L]. In Eig. we display the behaviour of the 
normalized and centered pair correlation functions gf{xi^Xj^t) (left) and gc{xi^Xj^t) 
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Figure 4. The normalized and centered pair correlation function gf^c{xi,Xj,t) 
for an initial collection of N = 10^ branching Brownian motions with diffusion 
coefficient D = 10 ~^ and birth-death rate A = 1 in a one-dimensional box of 
half-size L = 1. We take xj = 0 and plot gf^c(,Xi,Xj,t) with respect to Xi at 
successive times t = 1 (blue squares), t = 10 (red circles) and t = 40 (green stars). 
Symbols correspond to Monte Carlo simulations with 10^ ensembles, solid lines to 
exact solutions (Eqs. [TT] and |14| respectively). Left. For the case of a free system, 
gf{xi,Xj,t) initially develops a peak at Xi = xj = 0, which is the signature 
of particles undergoing spatial clustering. At later times, gf(xi,Xj^t) takes an 
asymptotic spatial shape, and is translated upwards by a spatially uniform term 
growing linearly in time. Right. For the case of a system with population 
control, gc{xi, Xj,t) initially develops again a peak at Xi = Xj = 0. Because 
of particle number conservation, an increased correlation about Xi = 0 implies 
negative correlations close to Xi = ±L. At l ater ti mes, gc{xi, Xj ,t) converges to 
an asymptotic spatial shape gT{xi,Xj) (Eq. |A.13|), displayed as a black dashed 
curve. 


(right) at successive times (their explicit expressions are reported in Appendix A). We 
set Xj = 0 and plot a function of —L < Xi < L. Solid curves represent 

the exact results given in Eqs. and respectively, at three increasing times t = 1, 
t = 10 and t = 40. Symbols represent Monte Carlo simulations performed with 10^ 
ensembles of 10^ neutrons. 

For the free system, the pair correlation function gf{xi^Xj^t) has three distinct 
regimes. Immediately after the initial time, gf{xi^Xj^t) displays a peak at short 


distances Xo 


which mirrors the effects of local fluctuations responsible for 


spatial clustering. The a mplitude of the peak is proportional to the ratio tdIte oc 
XL‘^/{ND) (see Eq. A.8), which precisely reflects the competition between migration 
and reproduction: the amplitude is larger for larger D and smaller A (for fixed L and 
A^), and vice-versa. The width of the peak, which is related to the correlation length of 
the system, is governed by diffusion, and is a growing function of D. For the limit case 
of non-diffusing particles (T) ^ 0), gf{xi^Xj^t) would display a delta-like behaviour at 
Xi = Xj, as expected on physical grounds: for long times, all the descendant particles 
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Figure 5. The average square distance between particles (r^) f,c{t) for fho 
one-dimensional model with N = 10^ initial neutrons, X = 1, D = 10“^ and 
L = 1. The blue solid curve corresponds to the free case: at long times, {r‘^)f(t) 
asymptotically converges to the ideal average square distance = (2/3) 

for a spatially uniform population, which is displayed as a blue dashed line. The 
red solid line corresponds to the case of population control: at long times, (r^}c(t) 
asymptotically converges to the value {r‘^)T given in Eq.which is displayed as 
a red dashed line. 


have died, except for a few point-like clusters composed of a very large number of 
individuals. For times shorter than the mixing time the amplitude of the peak 
grows due to births and deaths dominating over diffusion, whereas its width increases 
due to diffusion. When t >td^ the particles have explored the entire volume, and the 
tent-like shape of gf{xi^Xj^t) freezes into its asymptotic behaviour (see Eq. A.6). 

The total number of neutrons in the reactor also undergoes global fluctuations due 
to the absence of population control and to N being finite. These global fluctuations 
progressively lift upwards the shape of gf{ xi, by a spatially flat term that diverges 
linearly in time as ^ Xu 2 tlN (see E q.|A.8| ). Einally, for times larger than the extinction 
time te, gf{xi^Xj^t) > 1 (see Eq. A.7). This physically means that, no matter how 
dense the system is at time t = 0, global spatial fluctuations affect the whole volume 
with uniform (and increasing) intensity, and the neutrons are eventually doomed to 
extinction within a time t N/{Xu2) in the absence of population control (see Eig. 
b). 


The average square distance between particles for the free system is displayed in 
Eig.at time t = 0, the population is uniformly distributed and (r^)/(0) = {r‘^)id = 
(2/3)1/^. Immediately afterwards, {r‘^)f{t) starts to decrease due to spatial clustering. 
Eor times longer than global fluctuations dominate, and correlations range over 
the whole box. Then, {r‘^)f{t) increases and asymptotically saturates again to the 
ideal average square distance: this can be understood by observing that hf{xi^Xj^t) 
becomes spatially flat for t ^ te^ 

In the case of the system with population control, the pair correlation function 
gc{xi^Xj^t) has two distinct regimes. Immediately after the initial time, spatial 
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Figure 6. Monte Carlo simulation of the evolution of a collection of branching 
Brownian motions in a two-dimensional box, subject to population control. The 
particles are prepared at t = 0 on a uniform spatial distribution. In case a), 
the ratio between the migration area and the the specific square separation 
distance between neutrons is small, namely. A/{{r‘^)id/N) ~ 1.5, and clustering 
phenomena dominate over diffusion (however, since the total particle number is 
preserved, the population can not go to extinction). In case 6), the ratio between 
the migration area and the the specific square separation distance between 
neutrons is large, namely, A/{{r‘^)id/N) ~ 15, and spatial fiuctuations are much 
milder. 


clustering effects are again reffected in a peak at short distances Xi Xj for gc{xi^Xj^t). 
The amplitude and the width of the peak have the same behaviour as for the free 
system detailed above. However, due to the conservation of the number of particles, 
the positive correlations at the center of the box imply now negative correlations 
close to the boundaries Xi = ±L. For times shorter than the mixing time td, the 
amplitude of the peak grows and its width increases as in the previous case. Global 
spatial ffuctuations are intrinsically suppressed by N being fixed due to population 
control. For times larger than td, g^Xi^Xj^t) converges to an asymptotic spatial 
shape g^{xi,Xj) (see Eq. A. 13). In this regime, the amplitude of the pair correlation 
function is bounded by (see Eq. |A.14| ) 


\gc{xi,Xj,t)\ < 


(, 2 ) 


id 


NA 


( 8 ) 


where A = D/X is the characteristic migration area of the particles, i.e., the square 
distance explored by diffusion during a generation, and {r‘^)id/N is the specific square 
separation distance between neutrons corresponding to a uniform spatial distribution 
within the finite box [25]. In order for the fiuctuations to be small and prevent the 
emergence of spatial clustering, we must therefore have A':^ {r^)id/N, which occurs 
when the typical separation between particles is thoroughly explored within a single 
generation (see Eig. for a numerical illustration). Observe that the equilibrium 
condition for a reactor to be operated at the critical point does not depend on the 
total neutron population N: therefore, in a system with population control, spatial 
clustering can be quenched by simply imposing that N is sufficiently large, namely, 
N » {r^)id/A. 

The average square distance between particles for the system with population 
control is displayed in Eig. at time t = 0, the population is uniformly distributed 
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and we recover (r^)c(O) = {r‘^)id = (2/3)1/^. Immediately afterwards, {r‘^)c{t) starts to 
decrease due to the competition between diffusion and birth-death. For times longer 
than td, {r‘^)c{t) converges to the asymptotic value 


(r2)-= lim(r2),(t)=4F 

t^OO Arr) 


1 - 


I 2D 


tanh 


IXpL‘^ 

2D 


(9) 


which generalizes to confined geometries the findings for unbounded domains derived 
in m (see [Appendix A|) . 


4. The pair correlation function 

The function hf((Ki^^jD) for an exactly critical free system has been derived in [30] . 
and can be written as h/ = h^p + where (x^,Xj,t) = h^(xi,Xj,t) is the 
contribution from uncorrelated trajectories, and 

\u 2 N f dt' f dx'0(x^,x',t — t')0(xj,x',t — t')p(x',t') (10) 

Jo Jv 

is the contribution of the trajectories leading from the final positions at x^ and Xj 
at time t to the fission point x' at time t'. The coefficient 1^2 = ~ is 

the mean number of pairs created at each collision [36]. The term \i'2dt Np{yiD)d^ 
therefore represents the average number of ordered pairs created about x in the time 
interval (t,t + dt), which thus contribute to the correlation function about x^ and Xj 
at time t. Imposing a uniform spatial distribution Q = 1/V finally yields 

N(N-l) 

hf{yii,yij,t) = - ^ - -^\v2N J dt'Q{yii,yij,2t'). (11) 

The integral of the propagator appearing in Eq.[^is unbounded, so that at long times 
the amplitude of the correlations is expected to diverge. 

The pair correlation function he can be computed by closely following the 
arguments developed in uni. Actually, the reactor model described above is basically 
identical to that proposed in m, but for boundary conditions (neutrons evolve in 
a confined geometry, whereas in [T0| the viable space was unbounded) and initial 
conditions (in CO], all the particles were located at the same point at t = 0, whereas 
here the spatial distribution p of the individuals at t = 0 is arbitrary). Let us choose 
a pair of (distinct) neutrons located at x^ and Xj at time t. These neutrons may, or 
may not, have had a common ancestor (from a branching event) at a previous time 
0 < t' < t. Because of particle number conservation, the fraction of new particle pairs 
created in the time interval (t', t' -^dt) is Xpdt = Xi' 2 dt/{N — 1), obtained as the ratio of 
the new particle pairs created in the time interval, i.e., Xv 2 Ndtj 2 , to the total number 
of pairs N{N — l)/2 (TOj. The probability for a chosen pair of particles at time t not 
to have had a common ancestor is U{t) = so that the probability density for 

the ancestor to occur at time t' for a particle pair observed at t is 

Mt') = ( 12 ) 

The functional form of '0t(t') differs from that in [10], since the initial conditions are 
different. The function hc(x^,Xj,t) can be again written as he = hc^^ + hc^\ where 
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Xj,t) = U{t) hjjy^i^^j^t) is the contribution of neutrons having evolved freely 
with no common ancestors. The correlated contribution reads 


A^(A^ - 1) f dt' f dx.'Q{-Ki,x.',t - t')Q{-Kj,x.\t - t')'ipt{t')p{'x.\t').(13) 
Jo Jv 

where N{N — l)'ipt{t')dt' is the number of ordered particle pairs at time t that have a 
common ancestor in the time interval (t', T -\-dt'). The pair correlation function finally 
yields 


hc{^i,yij,t) = 


N{N-1) 




[ dt'e t/(xi,Xj,2i')(14) 
Jo 


when imposing Q = 1/V. The integral of the propagator appearing in Eq. is 
bounded thanks to the exponential term, and at long times the correlation function 
converges to an asymptotic shape. 

The specific shape of the correlation functions depends of the propagator, which 
can be computed once the dimension, geometry and boundary conditions of the 
problems have been assigned. In this respect, Eq. 14 generalizes the result by m 
in that it allows for arbitrary geometries and initial conditions. The average square 
distance can be then obtained by direct integration by following Eq. In the 
Appendix, we develop the calculations for one-dimensional domains, which can be 
easily generalized to arbitrary dimension. 


5. Conclusions 

In this work, we have analysed the statistical behaviour of fission chains in a prototype 
model of nuclear reactor. The neutron population evolving within a confined domain 
may display spatial clustering effects, due to the competition between births, deaths 
and diffusion. In particular, in the absence of population control, a reactor operated 
in stationary conditions at the critical point (corresponding to an exact equilibrium 
between the birth and death rates) would in principle lead to an almost sure extinction 
of the entire neutron population: this phenomenon goes under the name of critical 
catastrophe. Actually, such extinction is not observed in practice, because nuclear 
reactors are operated under strict population control policies (e.g., human intervention 
via control elements) that act a regularizing mechanisms at the global scale. 

A key ingredient for our analysis of the space-time behaviour of the fluctuations 
is the pair correlation function, which is proportional to the joint probability density 
of finding particle pairs at two spatial sites at a given time. We have in particular 
explicitly derived and compared the pair correlation function of a spatially confined 
neutron population that is left free to evolve to that of a neutron population that 
is kept under control by demanding that the total number of individuals is constant 
at any time. We have shown that in the former case the pair correlation function 
at the critical point diverges linearly in time, whereas in the latter converges to an 
asymptotic spatial shape. Eor a free system, the ultimate fate is therefore extinction, 
no matter how dense the neutron population is at the initial time. Eor the system with 
population control, spatial fluctuations can be tamed by acting on the total number of 
neutrons that are present in the reactor at equilibrium: when the ratio of the specific 
square separation distance between neutrons corresponding to an ideal uniform spatial 
distribution and the neutron migration area (i.e., the typical square distance travelled 
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by a particle by diffusion in the time span of a generation) is small, diffusion dominates 
over births and deaths, and the spatial clustering is quenched. This physically means 
that imposing a global feedback on the whole neutron population has a stabilizing 
effect also at the local scale of the spatial ffuctnations. 

The proposed reactor model retains the essential statistical aspects of actual 
systems, and as such it offers valuable insight on the behaviour of the spatial 
correlations. In our derivation, we have nonetheless introduced a number of simplifying 
hypotheses: for instance, we have assumed that the energy dependence of the physical 
parameters involved in neutron transport (such as the probability of fission and 
absorption, or the reaction rate) could be neglected; furthermore, we have modelled 
the boundaries of the reactor as refiecting, although a more realistic description would 
imply the possibility of particles leaking from the external surface. Finally, the effects 
of control rods in nuclear reactors are clearly localized in space, whereas in this 
work we have represented population control as a spatially homogeneous mechanism. 
Future work will be therefore aimed at investigating the impact of the introduced 
approximations. 
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Appendix A. Analysis of one-dimensional domains 


Consider a finite box of half-size L, i.e., a segment [—L^L] with V = 2L. At the 
boundaries x = TL, we impose refiecting (Neumann) conditions. The propagator for 
this system reads m 

1 1 

g{x, xo,t) = — + - '^ tpk{x)ipk{xo)e~“’‘\ (A.l) 

^ k=l 


where we have set 


and 


fk7r{L-x)\ 
cpk{x) =cos ( - ^ - 1 

.2 D 






(A.2) 


(A.3) 


We can identify the mixing time with td = (2/7r)^(L^/D). If we choose the uniform 
spatial distribution Q{xo) = 1/2L at time t = 0, the average density simply reads 


{n{xi,t)) = 


(A.4) 


As for the pair correlation function, the case of the free system is obtained by resorting 
to Eq. [E which yields 


hf{xi,Xj,t) = 


N{N - 1) Xiy2N 
(2L)2 “^(2L)2 


t + '^(pk{Xi)(pk{Xj) 


1 — e 


-2akt 


k=l 


Oik 


, (A.5) 
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For times t td, the exponential term in Eq. |A.5| vanishes, and the spatial shape 
of hf{xi^Xj^t) is frozen. In particular, the series appearing at the right-hand side is 
bounded, namely. 


's^^Tk{xi)ipk{xj) ^ 2 


k=l 


Oik 


3 D 


(A.6) 


-2k^-^ 




Then, for A" ^ 1 the amplitude of the normalized and centered pair correlation 
function asymptotically grows as 

^ {n{xi,t)){n{xj,t)) N ^ ^ 

Therefore, for times t ^ A/(Az/ 2 ), gf{xi^Xj^t) ^ 1, which allows identifying the 
extinction time te = N/ (Az^ 2 )- Observe that, for A ^ 1, gf{xi^Xj^t) can be expressed 
in terms of the two characteristic time scales, namely, 

t T ^ 1 — e 

gf{xi,Xj,t) = -^- (pkixi)ifik(xj) - 

The average square distance between particles 

2 , ^ J dxj J dxjjxj - Xjrhf{xi,Xj,t) 

^ f dxi f dxjhf{xi^Xj^t) 

can be obtained by integration. At time t = 0, 

(r"),(0) = 

At times t ^ te^ hf{xi^Xj^t) becomes spatially flat, and {r‘^)f{t) converges again to 
the ideal average square distance, namely, limt^oo(^^)/(^) = id' 

For the case of population control, from Eq. we get 

NiN - 11 Az/oA ^ 1 - e“0«fc+Ap)t 

hc{xi,Xj,t) = (2L)2 -,(A.ll) 


(A.8) 


(A.9) 


(A.IO) 


k=l 


Oik T ^ 


where Xp = A/(A — 1). Assuming that te ^ te, for times t ^ td the series appearing 
at the right-hand side is bounded, namely. 


E 

k=l 


iPk{xi)Tk{xj) 

Oik ^ ^ 


2ApL2 

D 


coth 


< 


2ApL2 

D 


- 1 


An 


(A.12) 


For A 1 the normalized and centered pair correlation function at long times 
converges to 


oo/ ^ \ — IR Tk\Xi)^k\Xj) 

9c[Xi,Xj) 2^ P + ^ • 

^ k=l 2r£; 

In particular, its amplitude is asymptotically bounded by 
, , ,, Xv22E 

\ge{Xi,Xj,t)\ < 

where the absolute value is taken because gc{xi,Xj,t) can be negative. 


(A.13) 


(A.14) 
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As for the average square distance, at time t = 0 we have again (r^)c(O) = {t‘^) 
as expected. The asymptotic behaviour of {r^)c{t) at times t ^ rjj can be computed 
exactly, and reads 




lim (r^)c(i) 



1 - 



(A.15) 


In the limit of extremely large populations, we have limjv-).oo(?’^)c = (2/3)i^ and we 
recover the ideal case. 
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